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We present a new multiphase-field theory for describing pattern formation in multi-domain and/or 
multi-component systems. The construction of the free energy functional and the dynamic equations 
is based on criteria that ensure mathematical and physical consistency. We first analyze previous 
multiphase-held theories, and identify their advantageous and disadvantageous features. On the 
basis of this analysis, we introduce a new way of constructing the free energy surface, and derive 
a generalized multiphase description for arbitrary number of phases (or domains). The presented 
approach retains the variational formalism; reduces (or extends) naturally to lower (or higher) num¬ 
ber of fields on the level of both the free energy functional and the dynamic equations; enables the 
use of arbitrary pairwise equilibrium interfacial properties; penalizes multiple junctions increasingly 
with the number of phases; ensures non-negative entropy production, and the convergence of the 
dynamic solutions to the equilibrium solutions; and avoids the appearance of spurious phases on 
binary interfaces. The new approach is tested for multi-component phase separation and grain 
coarsening. 


I. INTRODUCTION 

Despite recent advances in atomic scale continuum 
modeling of crystalline freezin^^ and efforts relying 
on the orientation field models®B2 the phase-field 
theoretical methods based on the multiphase-field 
(MPF) concept remain the method of choice, when 
addressing complex polycrystalline or multiphase/multi¬ 
component problems, such as multi-component phase 
separation or grain coarsening. A common feature of 
these models is that the individual physical phases / 
chemical components / solid grains are described by 
separate fields u(r, t) = [ui(r, £), ^2(1*, £), ..., i//v(r, £)]. 
A variety of this kind of models is available in the 
literature ranging from the early formulations by Chen 
and YanjJ^, SteinbacfP^, Steinbach and PezzolsPl, Chen 
and coworkersluny v j a later descendants by Nestler and 
coworkers Moelans and coworkers 21 24 ^ to more 

recent developments by Folch and PlappP^^, Kim et 
Takaki et a/P^, SteinbacfP^, Ofori-Opoku and 
Provatad^l, Cogswell and Carter 31 ^, Bollada et a/P^, by 
Emmerich and coworkers^, and Kim et alW^. These 
models differ in important details that improve the 
individual models in various respects relative to the 
others. It is, therefore, desirable to compare them from 
a theoretical viewpoint, and identify the possible advan¬ 
tages / disadvantages they have relative to each other, 
to see whether a more general formulation that unifies 
the advantageous features can indeed be constructed on 
the basis of the work done so far in this field. 


In attempting to develop a consistent description of 
interface driven multi-domain dynamics, we need to first 
identify the criteria the models have to satisfy. A few 
of such criteria have already been formulated along the 
development of the MPF models: 

(i) The multiphase-field descriptions view Ui(r,t) as 

the local and temporal volume/mass/mole fraction of 
the component/grain, prescribing thus = 1- 

(This work concentrates exclusively on these MPF mod¬ 
els; the multi order parameter theorieJ- 1 1 ! 14 * 15 * 21 ^ 24 * 30 *, 
which do not require this criterion, will be addressed else¬ 
where.) 

(ii) A further natural requirement is that the physical 
results should be independent of the labeling of the vari¬ 
ables. This condition is termed the ” principle of formal 
indistinguishability” of the fields. 

(iii) The solution of the dynamic equations should tend 
towards the equilibrium solution obtained from the re¬ 
spective Euler-Lagrange equations (ELEs) based on the 
free energy functional, where the equilibrium solution 
minimizes then the free energy of the system. 

(iv) As time evolves the free energy of the total system 
should decrease monotonically (second law of thermody¬ 
namics). 

(v) It is an evident requirement that the formulations 
for different numbers of phases or grains should be consis¬ 
tent with each other; i.e., it should be possible to recover 
the respective models from each other, when adding or 
removing a new phase/grain/orientation. It has been 
suggested recently that the usual variational approach to 
the MPF model does not satisfy this conditior l 26 * 32 *. 
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(vi) Another fairly general requirement, formulated by 
several authors, is that (a) the two-phase planar inter¬ 
faces should represent a (stable) equilibrium, and should 
be free of additional phases. This requirement can be 
extended to the dynamics (b) as follows: if a phase is 
not present, it should not appear deterministically (in 
the absence of fluctuations) at any time. This is of¬ 
ten called the condition of ”no spurious phase genera¬ 
tion”. The applicability of this requirement makes de¬ 
pends on the problem addressed: In the case of grain 
coarsening, for instance, an uncontrolled appearance of 
new grains / orientations at the grain boundaries needs 
to be avoided. (In other problems, however, metastable 
structures 35 37 or precipitates^ may appear at the solid- 
liquid interface, requi ring m odels that are able to de¬ 
scribe such phenomena 39 ^.) This criterion for the ab¬ 
sence of a third phase has been enforced different ways in 
different models. For example, Folch and Plapp^ have 
defined the free energy surface so that in the binary equi¬ 
librium the system stays in a two-phase subspace of the 
respective Gibbs simplex. Bollada et al 32 \ in turn, used 
the mobility matrix to force the system to avoid the for¬ 
mation of the third phase irrespectively of the free energy 
surface (and thus the equilibrium states of the system). 


II. CRITERIA OF PHYSICAL CONSISTENCY 

In this section, we present and discuss mathematical 
formulations of criteria (i) to (vii) identified above in de¬ 
tails. 


A. Free energy functional formalism 

In the multiphase approach, the following local con¬ 
straint [criterion (i)] applies for the variables: 

N 

yy*M) = i • (i) 

i=1 

As result of this constraint Eq. © is not an order 
parameter model, since normally different order parame¬ 
ters capturing various aspects of symmetry breaking are 
coupled to each other via physical laws, whereas here 
Eq. © prescribes a rather specific relationship: iq(r, t) 
represents the local fraction of the i —th phase, not 
identifiable as a quantity, whose magnitude is associated 
with the extent of symmetry breaking. 


Finally, a practical requirement: 

(vii) The model should allow the prescription of 
independent data for the interfacial properties and the 
kinetic coefficient of the individual phase pairs (including 
their possible anisotropy). 


While most of these criteria formulate natural/self- 
evident requirements, some of them were neglected when 
developing previous MPF models. 


In the present paper, we formulate an MPF approach 
that obeys all the criteria defined above. Herein, we 
address interface driven phenomena, in which the free 
energy density incorporates exclusively the ” interfacial” 
contributions, comprising the gradient energy and multi¬ 
well terms, whereas driving forces associated with tilting 
functions 26 are not considered. (An extension that in¬ 
cludes the latter will be outlined elsewhere.) The struc¬ 
ture of our paper is as follows. In Section II, we present 
the mathematical formulation of criteria (i) to (vii). In 
Section III, we first investigate which of these criteria are 
satisfied and which are not by the existing MPF models. 
We also point out which features of the individual ap¬ 
proaches can be adopted in developing a consistent the¬ 
ory. In Section IV we outline the generalized MPF formu¬ 
lation (henceforth abbreviated as XMPF) that satisfies 
criteria (i) to (vii). In Section V, we perform illustrative 
simulations using the XMPF approach to demonstrate 
the robustness of the theory. Section VI is devoted to a 
comparison with other models. Finally, in Section VII, 
we offer a few concluding remarks. 


The general interface contribution of a multiphase free 
energy functional is usually written as: 

F[u] = jdvtj Y^(V«i-V^)+^5(u)| , (2) 

where u(r, t ) is the vector of the variables, g(u) is the 
free energy density landscape, and A is a coefficient 
matrix of the general quadratic term for the gradients. 
For example, choosing A = I (where I) is the identity 
matrix yields a simple sum of the gradient square 
terms, A = I — e <g) e [where e = (1,1,..., 1)] results 
in a pure pairwise construction, while A^i = JV^?/ 2 , 
Aij^i = -UiUj corresponds to the anti-symmetrized 
(Landau-type) gradient term. The (generally) u- 
dependent coefficients e 2 and w can be related to the 
pairwise interfacial properties (the interfacial free energy 
and the interface thickness). 

The equilibrium solution can be found by solving the 
multiphase Euler-Lagrange equations: 

= A(r) , i = 1... N , (3) 

where 5F/5ui is the functional derivative of the free en¬ 
ergy functional with respect to iq( r), and A(r) is a La¬ 
grange multiplier emerging from the local constraint de¬ 
scribed by Eq. ©. Eliminating the Lagrange multiplier 
results in 


5F 

§Ui 


5F 

Su-j 


for = 1... TV , 


( 4 ) 
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thus offering the most general description of equilibrium. 


For non-conserved variables, the dynamic equations 
are written in the following general form: 


diLi _ j 

~dt~^ ij s^ 

3 = 1 J 


( 5 ) 


where the mobility matrix can be determined from the 
following conditions: 

1. The time derivatives sum up to zero [it follows from 
criterion (i)]: 


N 


E dui 

wr 1 = 0 . 

dt 


i= 1 


2. The variables are not labeled, i.e., none of them 
is distinguished on the basis of its index [criterion 

(ii)]- 

3. The solutions of the Euler-Lagrange equations must 
be stationary solutions of the dynamic equations [a 
requirement that follows from criterion (iii)]: 


<9u 

dt 



where u*(r) stands for a solution of Eq. Q. 

Applying condition 1 for Eq. 0 results in the general 
form 


dui ^ / or or \ 


N 


dt 


j=l 


SF SF 


Sm SlLo 


(6) 


where > 0 still can be arbitrary. Furthermore, using 
condition 2 yields 


N 




( 7 ) 




for j = 1... TV. Note that Eq. © and 0 resulted 
in a mobility matrix, whose elements sum up to zero in 
each row and column. Finally, condition 3 results in a 
symmetry condition, i.e. 


equations is to establish a free energy minimizing behav¬ 
ior [criterion (iv)]. Using Eq. (J5|, the condition for the 
time derivative of the total free energy reads as 



= -/ dV {(<5F u )L(( 5F u ) t } < 0 , 

where L is the mobility matrix, and we use the short no- 
tation SF U = (^-, ..., ^). Since the free energy 

must decrease in any volume, we can write 

(5F u )L(5F u ) t > 0 , (9) 

indicating that the mobility matrix must be positive 
semideftnite. A special, frequently made choice for the 
mobility matrix is the Lagrangian mobility = 1 /TV. 
Although this matrix is used quite widely, the dynamic 
equations are necessarily TV — dependent [i.e., criterion 
(v) is not satisfied], and it does not solve the problem of 
spurious phase appearance in non-equilibrium processes 
in general, as will be discussed later. 


B. Spurious phases 

Herein we address pattern coarsening phenomena, for 
which the requirement of ’no spurious phase appearance’ 
needs to be satisfied [criterion (vi)]. In general, this cri¬ 
terion means that any p-phase equilibrium solution (i.e., 
when exactly p fields are present) must be stable against 
the appearance of a new phase. Accordingly, assuming 
that q = TV — p of the TV phases (namely, ii, • • • ,i q ) 
are missing in an equilibrium solution u*(r), the condi¬ 
tion can be re-formulated as: 

5F = F[u*(r) + Su(r)] - F[u*(r)] > 0 

for any small perturbation for which Ylk =l ^ u k = 0 an d 
at least one of Sui 1 ( r), £tq 2 (r),..., Sui (r) is not equal to 
zero. In other words, leaving the p = TV — q dimensional 
subspace (together with keeping the local constraint, nat¬ 
urally) always has to result in higher energy. This condi¬ 
tion is satisfied for the equilibrium solutions representing 
minima of the free energy functional: Since 

F[u*(r)+<Su(r)] = F[u*] + — -Su T + ---+ ... , 


— ftji (8) 

(for the derivation, see Appendix A). Since 
K a — TV(TV — l)/2 mobilities can be 

chosen arbitrarily. 

Next we have to test the mobility matrix against the 
time dependence of the total free energy. The main rea¬ 
son of applying a linear approximation for the dynamic 


the second term on the right hand side vanishes 
for a solution of the Euler-Lagrange equation: 
(SF/Su) • £u = A(r) J2iLi $ u i = 0- Therefore, if B is 
positive definite, the equilibrium solution is a minimum. 
Consequently, if the binary planar interfaces represent 
minima of the multiphase functional, they are stable 
against the appearance of additional phases, which is 
a crucial requirement from the viewpoint of criterion (vi). 
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To fulfill criterion (vi), we also have to ensure the 
proper dynamic behavior of the system. The condition 
of ’no spurious phase appearance’ can be generalized for 
the non-equilibrium regime as follows: If a phase is not 
present in the system, it must not appear deterministi¬ 
cally, i.e. 


duj 

dt 


= 0 . 

Ui (r,t)=0 


( 10 ) 


Unfortunately, in an arbitrary , p-phase non-equilibrium 
state u(r), the condition SF = F[u(r) + 5u(r)] — 
F[u(r)] > 0 cannot be satisfied for all possible pertur¬ 
bations. Therefore, Eq. (10) cannot be guaranteed on 
the level of the free energy functional in general. Note, 
however, that the mobility matrix defines a ’conditional 
functional derivative’, which allows the system to leave 
the p-dimensional subspace only in particular directions 
[or, in other words, not any u(r) + Su(r) state is available 
from u(r) in the dynamics]. For special mobility matri¬ 
ces, like the Lagrangian matrix, one may find such a free 
energy functional that satisfies Eq. (10) 26 . Nevertheless, 
the free energy functional should not depend on the form 
of the mobility matrix in general. This problem is re¬ 
solved in a recent workP^l, in which the authors choose a 
mobility matrix having vanishing rows (and columns) for 
the fields not being present. Although this concept triv¬ 
ially results in Eq. (10), the application of such a matrix 
together with a particular free energy functional can be 
’dangerous’ in the sense that it may generate stationary 
solution from a non-equilibrium state, a possibility that 
has to be checked for all solutions obtained. 


III. ANALYSIS OF PREVIOUS MPF 
DESCRIPTIONS 

In this section, we analyze the most frequently used 
multiphase-field theories from the viewpoint of equilib¬ 
rium solutions. We check whether the trivial extension 
of the planar interface emerging from the binary reduc¬ 
tion of the free energy functional is a solution of the 
multiphase problem too. Summarizing, we require that 
the equilibrium solution for the pure binary planar inter¬ 
face in the multiphase-field problem (N < 3) coincides 
with the equilibrium solution for the binary planar in¬ 
terface of the binary (N = 2) problem [criterion (iv)]. 
The methodology for testing this feature consists of the 
following steps: 

1. Take the free energy functional in the two-phase 
limit; 

2. Solve the respective Euler-Lagrange equation of the 
two-phase problem for planar interface geometry, 
resulting in u(x)\ 

3. Make a natural multiphase extension of u{x) via 
adding zero additional fields as needed for the mul¬ 
tiphase case, i.e. Ui(x) := u(x), Uj^i(x) := 1 —u(x), 
and Uk^ij(x) = 0, where 1 < i, j, k < V; 


4. Plug the extended solution into the Euler-Lagrange 
equations of the multiphase problem described by 
Eq. ([3]), and check whether it satisfies them. 

The test can be simplified in case of free energy func¬ 
tionals constructed exclusively from pairwise contribu¬ 
tions, i.e. 



where J2i<j = ZliLi 1 Z^L*+i • f(u,v) is called genera- 
tor. The functional derivatives read as: 

^7 = Yl S K u i> u i) > (U) 

1 


where 


Sf(ui,Uj ) 


df v df 

dui Vdui 


Assuming that 5f(ui, 0) = 0, and plugging in the ex¬ 
tended planar interface solution into the Euler-Lagrange 
equations of the multiphase problem yield 

= Sf[u(x), 1 - u(x)} = X(x) 

OUi 

J— = <V[! - u(x), u(x)} = X(x) 
duj 

4— = */[0, if(z)] + <5/[0,1 - u(x)) = X(x) , 
du k 

where Ui(x) = u(x), Uj^i(x) = 1 —u(x), and Uk^ij(x) = 
0, and u(x) represents the planar binary solution. From 
the equations above follows that 

5f(u, 1 - u) = <5/(1 - u, u) = <5/(0, u) + <5/(0,1 -u) (12) 

must apply. If Eq. ( |12| ) does not apply, the binary planar 
interface is not an equilibrium solution of the multiphase 
problem. 


A. Steinbach et al. 

In 1996, Steinbach et al. 12 proposed the following free 
energy functional for multiphase systems, which serves 
as the basis for the worldwide used phase-field software 
MICRES^l 

F = J dV {/i„ t f(u, Vu) + f di (u)} , (13) 

where Y^=i u i{ Y ^) = whereas the term /df(u) is 
responsible for the thermodynamic driving force (i.e. 
free energy difference between the bulk phases) and 
/intf( u , Vu) denotes the interface energy consisting of a 
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gradient [/ gr (u, Vu)] and a multi-well [/ mw (u)] contribu¬ 
tion: 


/intf (U, Vu) = / gr (u, Vu) + fmw (U) , (14) 

where the terms are given in the following specific forms: 

e 2 - 

/gr(u, Vu) = 22 f ( M *VMj - UjVUi) 2 
i<j 

/mw(u) = Yj ~2~( u i u j) 2 • (15) 

i<j 

The functional naturally reduces to the standard binary 
form 


F = 



Wij_ 

2 


N 1 - «)] 



which are variational, therefore, the planar binary inter¬ 
faces are not a stationary solutions of these. To avoid 
the problem, the authors used a ’’binary approximation” 

, in which all terms of k ^ i,j indices are 


of 6F_ _ SF_ 
Ui SUi Su-j 

neglected 12 : 


SF dF 

Sui Suj 


Wi 


UiUj{uj — ui) + eUui^Uj — Uj\7 2 Ui ) 


Although the planar binary interfaces represent station¬ 
ary solutions of the resulting non-variational dynamics, 
Eq. © does not apply, therefore, the dynamics is not en¬ 
ergy minimizing in principle , as it will be demonstrated 
later. 


B. Steinbach-Pezzolla 


The ID Euler-Lagrange equation reads as: dF/du = 
Wiju{ 1 — u)( 1 — 2 u) — ejjd 2 u = 0, thus resulting in the 
usual 

_ 1 + tanh[a;/(2<%)] 
u{x) — ^ 


In the Steinbach-Pezzolla formalism (published first in 
1999^, and adopted in various workJ^ 29 31 including the 
OpenPhase softward^), the interface contribution to the 
free energy is given by the following simplified form of 
Eq. (pb: 


planar interface solution, where d 2 - — e 2 - / Wij . As a first 
step, we investigate whether the extension of this solution 
minimizes the multiphase problem. Since fi n tf{ u, Vu) is 
a pure pairwise construction, it is enough to take the 
generator, which reads as: 

f(Ui,Uj) = -fimVUj - UjVUi) 2 + y- u 2 u 2 , 
yielding 


fin £ = > ( 2 °) 

i<j 1 ' 2 

where the gradient term —Viq • Vuj is the linear ap¬ 
proximation of (UiVuj — UjVui ) 2 , <Jij the interfacial free 
energy of the equilibrium (i,j) planar interface, while the 
local term \uf \uj\ is responsible for the finite interface 
width given by g. Reducing the model to the N = 2 case 
yields 


df (ui) Uj) — WijUiUj + Cij ^2(uid x Uj Ujd x Ui)d x Uj+ 

+ ( Uid 2 x Uj — Ujd 2 Ui)uj] . 

(16) 


Note that df(u, 0) = 0. Substituting Ui(x) = u(x), 
Uj^i(x) = 1 — u(x) and Uk^ij(x) = 0 into Eq. (16), 
yields 


F- = Sf(u, 1 - u) = y-sech 4 [a;/(2(5 i j) 

y- = <5/(1 ~u,u) = y-sech 4 [;r/(2 V,) 
5F 

-— = (5/(0, u) + (5/(0,1 - u) = 0 . 
duk 


(17) 

(18) 
(19) 


These equations clearly show that the planar binary 
interfaces are not equilibrium solutions of the multiphase 
problem. 

The dynamic equations read as: 

_dui _ y /SF SF\ 

dt y* u \ S'U j 5uj) ’ 


SP _ ^. LJ 


4 Gi 


\^Ju) 2 + \u\ \l-u\ 


( 21 ) 


where we used u(x) := Ui(x ), Uj(x) := 1 — u(x), and 
Uk^ij(x) = 0. The ID Euler-Lagrange equation reads as 

f2g 2 \ 

sign (u) |1 — u\ — sign(l — u ) \u\ = ( —y I d 2 u{x) , 


from which 


l(x) = 


1 + sin[(7r/?7) • x] 


emerges for the planar binary interface. The generator 
function reads as 


^ 4(7 i ' f TJ 2 

f (Ui, Uj) = I \Ui\ \Uj\ • Vtq 

g \ 7T 


■3 > 


therefore, 


5f(ui,Uj) = ^sign(uj) \uj \ + y V 2 Uj 
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yielding 

- 2 

Sf(u, 1 - U) = -CTij 

- 2 

Sf(l-u,u) = -dij 

<5/(0, u) + <5/(0,1 - u) = - sin i-x] (o-fej - o-fci) 

V \V / 


for the extension Ui{pc) = iz(x), Uj =7 ci(x) = 1 — 'u(x) and 
Uk^ij(x) = 0 , showing that the planar binary interfaces 
do not minimize the free energy functional for non-zero 
interfacial free energies. 


The dynamic equations of Steinbach and Pezzolla 13 
also read as: 


duj 

dt 




'6F_ 

Sui 


SF ' 
S Uj _ 


with Kij — Kji > 0 constants, therefore, the mobility 
matrix satisfies conditions 1-3 of Section II.C. Consider¬ 
ing that the planar binary interfaces represent no equi¬ 
librium solution, they are not stationary. The prob¬ 
lem is resolved again by replacing the gradient term 
—'Vui-'Vuj by ( Ui'Vuj — Uj\7ui ) 2 , and using the ’binary 
approximation' 29 : 


SF 5F 

Sui Suj 


sign(rq)|rtj| - sign(^)|^| 


— ( UjAui — UiAuj ) . 


( 22 ) 


functional. In addition, the absence of the absolute value 
function terminates the bulk rq = 1 equilibrium solution 
too. 


C. Nestler-Wheeler 


Another descendant of the original Steinbach et 
al. models is the general Nestler-Wheeler type 
formalismP^® 


/NW _ 
J intf 


i<j 


ij 


2 {uiVuj - UjVui ) 2 + ^-(K 


\ u i\Y 


{23) 

where p = 1 or 2. For p = '/HDEH] gq (23) recov¬ 
ers E qs. (|T3| )-(|l5| (Steinbach et al.), whereas in case of 
p = ilizii9Enr ^ reduces to Eq. ( 21 ) for N = 2 with the 
solution 


, , 1 + sin(x/<5jj) 

«(*) = -g- • 

For p = 1, the derivative of the generator function reads 
as 


5f(ui,Uj ) =(wij/2) sign(ui)\uj\ + 

A ^ij [2 (uid x Uj UjD x Ui^O x Uj-\- 
+ (uid x Uj - Ujd x Ui)uj\ , 

which, in the case of Ui(x) := u(x), Uj^i(x) := 1 — u(x), 
and Uk^ij(x) := 0 , yields 


Although the resulting non-variational dynamics stabi¬ 
lizes the extension of the planar interface solution u(x) = 
{ 1 + sin[( 77 / 7 r)x]}, unfortunately it does not minimize the 
free energy functional. A remarkable improvement of the 
Steinbach-Pezzola model has been put forward by Kim 
et all 27 . Introducing step functions that are Si = 1 for 
Ui > 0 and Si = 0 otherwise, the mobility matrix has 
been assumed to have the form shown below 


duj 

dt 


E 

j#* 


^ij SiSj 


Su* 


AF’ 

Suj_ 


This change leads to an important step ahead: it 
retains the variational formalism, while stabilizing the 
flat interface. Nevertheless, it is not yet a solution 
of the Euler-Lagrange equation of the multiphase 
problem, therefore, this mobility matrix is ’dangerous’ 
in the sense that it generates stationary solution from 
a non-equilibrium state. This approach has recently 
been applied for an asymmetric case taking the grain 
boundary energies from a databasd^. 


Finally, we mention that the derivations presented 
above can be trivially repeated for using UiUj instead 
of \uj\ \uj\ in the free energy functional described by Eq. 
( 20 ), resulting in the same qualitative results, i.e. the 
planar binary interfaces do not minimize the multiphase 


A 3 ? if ■ 

5f(u, 1 - u) = -22- cos 2 {x/Sij) 

A 3? n •• 

<5/(1 -u,u) = -^p-cos 2 {x/Sij) 

<5/(0, u) + <5/(0 ,1 — u) = 0 , 


showing that the binary planar interfaces do not min¬ 
imize the free energy functional again. This, together 
with the fact that the mobility matrix was chosen to be 
Lagrangian, means that we have here the same problem 
as in the case of the Steinbach-Pezzolla formalism, which 
can be resolved by using Eq. (22), i.e., by adopting 
non-variational dynamics. 


In a recent variant of the p = 1 Nestler-Wheeler model 
by Ankit et the multi-obstacle free energy landscape 
contains a triplet term of the form: 

fz = E ^ijk\ui\\ui\\u k \ , (24) 

i<j<k 

where the triple sum runs for all different (i, j, k) triplets, 
and the authors use 7 ^ to control the appearance of 
the third phase at the binary interfaces. We note, 
however, that Eq. ( [24] ) has no effect on the existence 
of the planar interface solution, since the derivative 
dfs/dui = sign(rq) Xq</c \ u i\ \ u j\ vanishes for binary 
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planar interfaces. In other words, fy is not suitable 
for generating equilibrium planar interface solutions. 
Nevertheless, choosing 7^ —»• oc results in two-phase 
interfaces free of additional fields, but for any finite 7 ^, 
additional fields are always present at the interfaces 
mathematically. 

We note that, in models relying on the multi-obstacle 
potential, /i nt f oc is often prescribed out of the 
physical regime to prevent the evolution of the fields 
into the ” unphysical” states iq < 0 and rq > 1. This is 
another way to stabilize the (otherwise non-equilibrium) 
two-phase interfaces. We recall furthermore that there is 
physical interpretation for iq < 0 and iq > 1. Compari¬ 
son of the phase-field models to the Ginzburg-Landau 
model and/or to amplitude equations emerging from 
classical density functional theories^ implies that tq < 0 
can be simply associated with a negative amplitude of 
the first reciprocal lattice vector set in the crystal, which 
is a real perturbation of the liquid state for cubic crystal 
structures. Similarly, iq > 1 is nothing more than an 
amplitude larger than the equilibrium crystal amplitude. 


D. Folch-Plapp 


The term ” Lagrange multiplier formalism” originates 
from Folch and PlappPn In their multiphase description, 
the free energy functional is based on several theoreti¬ 
cal considerations including binary equilibrium solutions 
and the condition of no spurious phase generation. For 
three phases the interface contribution of the free energy 
functional reads as: 

/iSf = j E ( Vu *) 2 + jfrw(u) , (25) 


where the ” triple-well” free energy density is frw( u) = 
-where g(iii) = [u»(l - Ui)] 2 . First, we an- 
alyze the binary planar interfaces. For N = 2 the free 
energy functional reduces to 

fintf = £ 2 (V «) 2 +wg{u) . 
generating the usual ID Euler-Lagrange equation 

= wg'(u ) — e 2 d 2 u = 0 . (26) 


The equilibrium planar interface solution is then u(x) = 
{1 + tanh[x/(25)]}/2 with S 2 = e 2 /w. Since Eq. (25) is 
not a pairwise construction, the generator function tech¬ 
nique does not apply here. The general Euler-Lagrange 
equations read as: 


SF_ 

Sui 


\ [w g'(Ui) - € 2 d 2 x Ui] = A(r) 


(27) 


Comparing Eqs. (26) and (27) results in 2 important 
properties of the model: 


SF 5F 
Sm Su 


and 


Sui 



indicating that the equilibrium planar binary interfaces 
are solutions of the multiphase problem with X(x) = 0 . 
The proposed dynamics uses the Lagrangian mobility 
matrix , therefore, the Folch-Plapp model is the first 
model, which passes the binary interface criterion. 


Next, we discuss the appearance of spurious phases. 
The general condition reads as u k (r,t) = 0, i.e. if phase 
k is not present at t = 0 , it must not appear at any time. 
In the Folch-Plapp model, the time evolution of a phase, 
which is apparently not present reads as 


where k 7 ^ z, j and i 7 ^ j. Using Eq. ( |27| ) it is obvious 
that SF/Sui + SF/Suj s 0 for iq(r, t) + Uj(r,t) = 1, 
therefore, no spurious phases appear. 



du k f 5F 5F 

dt Uk= o ^ 5u i 


The authors have also worked out an asymmetric 
version of the model, in which different binary interfacial 
free energies can be used, which also satisfies the 
basic criteria of physical consistency together with no 
spurious phase generation. Despite its advantageous 
features compared to former approaches, there are a few 
weaknesses of the model: (a) to avoid the appearance of 
spurious phases, a free energy functional is used whose 
form depends on the particular choice of the mobility 
matrix, which is clearly not physical, (b) no N > 3 
generalization of the model is available. Furthermore, 
(c) as the authors suggested, the necessary minimum 
exponent in the multi-well term might be proportional 
to the number of phases, making the model practically 
useless when a large number (e.g., thousands) of dif¬ 
ferently oriented dendrites or grains have to be simulated. 

At this stage, it is clear that the Folch-Plapp model is 
a large step towards constructing a consistent multiphase 
description, since it satisfies almost all the criteria of 
physical consistency. Unfortunately, the price is high: it 
is not immediately clear how to generalize this model to 
N > 3, together with avoiding the appearance of spurious 
phases via introducing special terms in the free energy 
functional, which follow from the form of the dynamic 
equations. 


E. Bollada-Jimack-Mullis 

In the previous sections it has been demonstrated that 
the condition of no spurious phase generation works at 
two levels. First, if a solution of the Euler-Lagrange equa¬ 
tion represents minimum of the free energy functional, 








it is stable against small perturbations (assuming vari¬ 
ational dynamics with a positive semi-definite mobility 
matrix). In addition, to avoid the appearance of spu¬ 
rious phase outside of equilibrium may necessitate the 
adjustment of the free energy functional. In a recent 
work 32 , however, Bollada et al. avoided the problem by 
introducing a field-dependent mobility matrix 

N 

La = h(ui,iLj ) and Lij = —h(ui,Uj) for j , 

j/* 

(28) 

where 


h(ui,Uj) = 


1 — Uj 


1 — Un 


Apparently, Eq. (28) satisfies the conditions of Section 
II.C; i.e., the dynamics ensures non-negative entropy 
production inside the simplex (i.e., when all iq £ [0,1], 
all h(ui,uj ) >0). In addition, the spurious phase 
generation is excluded, since for Uk(r,t) = 0 the k th row 
(and column) of the mobility matrix vanishes, yielding 
(duk/dt)\ Uk ( r ,t)=o = 0- Note that this is achieved with¬ 
out revising the free energy functional. Moreover, the 
description is N- independent, since the mobility matrix 
consistently reduces to the N — 1 case. Despite the 
significant improvement, one has to be careful with the 
Bollada-Jimack-Mullis mobility matrices, since they can 
be dangerous with respect to the free energy functional 
in the sense that the mobility matrix may generate sta¬ 
tionary solution of the dynamics from a non-equilibrium 
state (as it happens in case of a multi-obstacle potential). 
In addition, outside of the simplex (i.e., for iq < 0), the 
mobility matrix is indefinite, thus violating criterion (iv). 


The results of the above analysis of the MPF theories 
are summarized in Table I. The following has been found: 

(A) The criteria for the sum of the fields and for label¬ 
ing are satisfied by all the models investigated. 

(B) The lack of equilibrium planar two-phase inter¬ 
faces (in the 7V-phase problem) can be resolved by em¬ 
ploying: (a) non-variational dynamics (the models by 
Steinbach et all 12 } Steinbach and Pezzol#^, Nestler and 
WheelePl®), (b) a degenerate mobility matrix (the 
models by Kim et a/. 2 ^, and by Bollada et a/. 32 ^), or by 
(cipenalizing the triplet term (the model by Ankit et 
aZP^j). We identified the following problems associated 
with methods (a)-(c): the solution does not converge 
to the equilibrium solution; furthermore, in case (c) the 
third phase is unavoidably present (even if in a small 
amount). 

(C) When the equilibrium conditions are satisfied (as 
in the model by Folch and PlappP2), we obtain an TV- 
dependent approach without the possibility of prescrib¬ 
ing freely the pairwise interfacial data. 

Considering these, we conclude that none of the MPF 
models investigated here satisfy all the consistency crite¬ 
ria specified. We stress furthermore that the introduction 


TABLE I. Properties of different multiphase-field models from 
the viewpoint of criteria defined in this work. 


model \ criterion 

i 

ii 

iii 

iv 

v 

vi( a ) 

vi(b) 

vii 

Steinbach et al. 

X 

X 





X 

X 

Steinbach & Pezzola 

X 

X 





X 

X 

Nestler & Wheeler 

X 

X 





X 

X 

Kim et al. 

X 

X 


X 

X 


X 

X 

Bollada et al. 

X 

X 


X 

X 


X 

X 

Ankit et al. 

X 

X 

X 

X 




X 

Folch & Plapp 

X 

X 

X 

X 


X 

X 



of additional thermodynamical driving force via an ap¬ 
propriate tilting function (as needed for describing poly- 
crystalline solidification) would influence neither the va¬ 
lidity of these criteria, nor the outcome of this analysis. 


IV. CONSISTENT MULTIPHASE FORMALISM 
A. General framework 

Herein, we derive a multiphase description that satis¬ 
fies criteria (i) to (vii). It is useful to start with the con¬ 
dition of formal reducibility [criterion (iv)]. First, setting 
?i]v(r, t) =0 in the V-phase free energy functional 
should result in the N — 1 phase functional, F N ~ 1 . The 
same should apply for the dynamic equations 

_du i _ sr fSF SF\ 

as well; i.e., -u^ = -5F^ | WiV ( r ,t)=o should reduce 

to — • SFu N_1 \ and un = 0. Note 

that the latter satisfies the condition of no spurious phase 
generation [criterion (ii)]), since = 0 = 0. HereL^ 
is a general, symmetric, positive semidefinite V-phase 
mobility matrix, i.e. = —Kij for z ^ j, while La = 
Yhjz /,1 where i, j = 1... V, and > 0. Since 

the (modified) Bollada-Jimack-Mullis matrix defined by 
Eq. (J28| satisfies the condition of formal reducibility, we 
choose this mobility matrix, namely, 


Ui 


Uj 

1 Ui 


1 -Uj 


where constant positive prefactors accounting for the 
mobility of the planar i, j interface are also incorporated. 
Furthermore, we prescribe the condition of reduction 
also for the functional derivatives; i.e., the first N — 1 
components of 5F^ should reduce to SF^ N ^ for 
t) = 0. Since the N th row (and column) of the 
reduced Bollada-Jimack-Mullis matrix is 0, it always 
results in un = 0, therefore, (SF N /Su]\f)\u N (r,t)=o can 
be arbitrary. Since the n-phase Bollada-Jimack-Mullis 
matrix is positive semidefinite on an n-phase state (i.e. 
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when none of the n components vanishes) with multiplic¬ 
ity 1 for the eigenvalue 0 (it can be proven numerically), 
the n-phase stationary solutions coincide with the n 
phase equilibrium states. Moreover, since the dynamic 
equations reduce naturally to the (AT — l)-phase case, 
the stationary solutions of the reduced dynamics include 
the natural AT phase extensions of the (.AT — l)-phase 
equilibrium solutions (where none of the N — 1 phases is 
missing). Therefore, the natural AT-phase extensions of 
all (AT — l)-phase equilibrium solutions emerging from 
p( N -i) represent extrema of F N . If this is true, 

the Bollada-Jimack-Mullis matrix is not dangerous with 
respect to the free energy functional, since all stationary 
states of the dynamics represent equilibrium. Since 
the condition must apply for arbitrary AT, the general 
condition for the free energy functional reads as follows: 


The (p + q)-phase trivial extensions of all p-phase 
equilibrium solutions (where all the p phases are 
non-vanishing) emerging from the p-phase free en¬ 
ergy functional F( p ) must represent extrema of the 
(p + q) -phase free energy functional F(p+ q ) too, for any 
p > 0, q > 0. 


For practical reasons, we introduce the following con¬ 
dition: for a field, which is not present, the functional 
derivative vanishes, i.e. 

5F_ 

Sui 



If the Lagrange multiplier also vanishes [A(r) = 0] for 
all p —phase equilibrium solutions of F ^ (where p > 1 
arbitrary), while the free energy functional (and the func¬ 
tional derivative) reduces naturally, all trivial extensions 
of all p-phase equilibrium solutions remain equilibrium 
solutions of the AT-phase free energy functional. Conse¬ 
quently, here the Bollada-Jimack-Mullis matrix does not 
stabilize non-equilibrium solutions, while preventing the 
appearance of spurious phases. Note, that this is obvi¬ 
ously not true for the models of Steinbach et a/., Stein- 
bach and Pezzola, Nestler and Wheeler, Ankit et a/., and 
for the model potential used in the work of Bollada, Ji- 
mack, and Mullis. Although there the free energy func¬ 
tionals and the functional derivatives reduce naturally, 
and (SF/5ui)\ Ui= o = 0 also applies, the planar two-phase 
interface solution generates different Lagrange multipli¬ 
ers for the vanishing and non-vanishing fields in the com¬ 
plete (AT-phase) Euler-Lagrange problem, indicating that 
the natural extensions of the planar two-phase interfaces 
do not represent equilibrium. Yet, the application of the 
Bollada-Jimack-Mullis mobility matrix transforms them 
into stationary solutions, showing that in these cases the 
application of the Bollada-Jimack-Mullis matrix is ” dan¬ 
gerous” . 


B. Free energy functional 

1. Symmetric system 


The main question is, how one should construct an 
interface term that satisfies the conditions given above. 
First, we consider the symmetric case, where all interface 
thicknesses and interfacial free ene rgies are equal. Fol¬ 
lowing Chen and co-workers 11 ^3E 15 and Moelans and co¬ 
workers 21 the interface term of the free energy func¬ 
tional is constructed as follows 

2 N 

/intf = J yyv«i) 2 +wg( u) , (29) 

i= 1 

where we use the following new Ansatz for the multiphase 
barrier function: 


p(u) 




2 2 
u iUj 


(30) 


The functional derivative reads as 


SF 


W [Ui( U 2 — Ui)\ — 6 2 \7 2 Ui , 


(31) 


which vanishes for iq = 0. The binary planar interface so¬ 
lution is u{x) = {l + tanh[x/(2J)]}/2 (where S 2 = e 2 /w), 
for which the multiphase Euler-Lagrange equations re¬ 
duce to 

= wu(l — u)(l — 2 u) — e 2 d 2 u = 0 
Sui Suj v ' v ' * 

Here we used the trivial extension rq := u(x), 

Uj^Lj := 1 — u(x) and u^ij = 0. Since the free 
energy functional and the functional derivatives reduce 
naturally, and SF/Sui vanishes for iq =? 0, together with 
the fact that any trivial extension of the planar interface 
solution represents equilibrium, it is reasonable to assume 
that the Bollada-Jimack-Mullis matrix is not dangerous 
considering at least the planar interfaces, i.e., it does 
not stabilize a non-equilibrium planar interface, since all 
planar interfaces represent equilibrium. Naturally, the 
same investigation should be repeated for all n-phase 
equilibrium solutions of F^ for any positive n, however, 
this kind of study is out of the scope of the present paper. 


It is important to mention, that Eq. (30) shows a very 
practical feature 


g({l/N,l/N,...,l/N}) = 


12 


1 - 


AT 2 


(32) 


i.e., the higher-order junctions are energetically increas¬ 
ingly less favorable. Note that this is not true for Eq. 
(15). The tendency of increasing free energy is also en¬ 


sured by the multi-well term defined in Eq. (20), however 
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there, as we have shown previously, the binary planar in¬ 
terfaces do not minimize the free energy functional in 
the general 7V-phase case. It is worth noting that Eq. 
(32) contradicts Folch and PlappPS, who expect that the 


polynomial degree of the g( u) function that penalizes the 
high-order multiple junctions would increase with N. Eq. 
(32) shows that the double-obstacle function is not the 


only one that realizes this tendency, furthermore, here 
[see Eq. (30)] the planar binary interface represents an 
equilibrium solution of the multiphase problem. 


Note that 

duk 


u k =0 


_ dx 

dUi n 


i~\~ u j —1 


= 0 , where x = 


e 2 ,re, therefore, (5F / 5uk ) Uk =o = 0 , and the second line 


of Eq. (36) also vanishes for iq = u, Uj = 1 — u , = 

0 , therefore, the planar binary interfaces are equilibrium 
solutions of the multiphase problem for arbitrary pairwise 
e 2 - and Wij fitted to the interfacial free energy 7 ^ and 
interface thickness Sij as follows: 


efj = 3 (S^ • 7 ij) and 


Wi 


— 3 ( 7 ij/Sij) 


2. Asymmetric extension 

Following Moelans- 2 , the asymmetric extension of 
Eq. ( [29] ) can be obtained by employing the Kazaryan- 
polynomialJ^ 


2 (u) := 




Ei 


w( u) := 


E- 


;3 * 3 

■i,s Wjju^uj 


E,, 

The free energy density then reads as 


o / \ N 

e 2 (u) 


(33) 

(34) 


/intf = —y](VMi) 2 + w(u) g( u) , (35) 


3. Introducing anisotropy 

In various practically important cases, including den¬ 
dritic solidification and grain coarsening, the interfa¬ 
cial free energy between two phases displays anisotropy, 
which can be formulated mathematically as: 

6 ij e ij I 1 + CLij • rjij(Uij)] , 

where aij is the amplitude (strength) of the anisotropy, 

Viq — Vuj 


| Wui — Wuj 


i= 1 


is a unit vector characterizing the (i,j) binary interface, 
while rjij(iiij) reflects the crystal symmetry. This exten¬ 
sion modifies Eq. (33) and the functional derivative as 
follows 


where g( u) is defined by Eq. (30). Since the Kazaryan 
polynomial is ” quasi-constant” (the nominator and the 
denominator contain the same terms with different co¬ 
efficients), it is reasonable to assume that this modifica¬ 
tion does not change the structure of extrema of g( u). 
Although it will be demonstrated for asymmetric N = 4 
and N = 5 systems, a strict mathematical derivation for 
arbitrary number of phases is out of the scope of the 
present paper. 

The binary reduction of Eq. ([35]) reads as 


/intf = 4(V«) 2 + Wij[u( 1 - It)] 2 , 

generating the equilibrium planar binary interface with 
<5 2 = e^/wij. The general functional derivative is 


5F dg dw 2 2 

7 — = w -- \- — g{ u) - e V Ui 

OUj 


dui 
N r 

E 

j=i 


dui 

1 de 2 „ de 2 „ 

\7Un — — - \/Ui 


2 duj 


8Un 


■ \/Ug 


(36) 


SF , ,dw 
5m ~ 9 ^ dui 


■ w(u) 


dg 

dui 


de 2 

d Ui 


1 N 

^E(v«i) 2 


3 = 1 


- V- 


de 2 

d\7ui 


N 

E(^ 

3 =1 


+ e 2 • Vui 


Here the extra term reads as 


^ 2 Em 


,2 n ,2 


d\7ui 


Em u 2 k u 2 


where 


det 


Mr- oc u 2 , which 


Uj = 1 


u IS 


dx/ui ^ therefore, 

means (SF/Suk)\ Uk =o = 0- In addition, for u, 

de 2 _ de % therefore 5F — i w heri 

dVuij dVu i}j > ^nereiore, Smj ~ 2 5u ’ wner 

the functional derivative in the reduced model, there¬ 
fore, the equilibrium binary interfaces emerging from 
SF/Su = 0 are stationary solutions of the multiphase 
problem. 


where 




dw 

M = { i] 


Td^ji^ij - wyuf 


Em u 2 u 2 


V. TESTING THE XMPF MODEL 

In this section, we review whether the proposed model 
satisfies indeed criteria (i) to (vii). Several of these crite¬ 
ria are satisfied owing to the specific formulation of our 
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model [these are (i), (iii) - (vi), and (vii)] as summa¬ 
rized in sub-section A. Fulfillment of practical criterion 

(ii), however, needs further investigation, which is under¬ 
taken in sub-section B. Finally, illustrative simulations 
are presented for grain coarsening in sub-section C. 


A. Consistency criteria satisfied 

It can be shown that the proposed model satisfies the 
following criteria: 

(i) J2iLi u i( r ^) = 1 ( see Section IV.A); 

(ii) Since the mobility matrix is symmetric, the 
physical results are invariant to formally exchanging 
pairs of field indices, i O j, i.e. the variables are not 
labeled (see Section IV.A); 

(iii) Under appropriate boundary conditions, any 
trivial multiphase extension of the equilibrium binary 
solution represents equilibrium solution of the multi¬ 
phase problem, which is then a stationary solution of the 
dynamic equations towards which the time dependent 
solution evolves (see Section IV.D); 

(iv) Since the mobility matrix is positive semidefinite, 
non-negative entropy production is ensured for both 
conserved and non-conserved dynamics (see Sections 
IV.A and VLB); 

(v) Reduction/extension of the TV-held theory to 
TV — 1 or TV + 1 fields is trivial on the level of both the 
free energy functional (and the functional derivative) 
and the mobility matrix (see Sections IV.A and IV.D); 

(vi) No additional phases appear at the equilibrium 
binary interfaces (see Sections IV.B, IV.C, IV.D and 
VI.A), and the dynamic spurious phase generation is 
also excluded (see Sections IV.A and V.B); 

(vii) Freedom for choosing independent inter facial (e^- 
and Wij ) and kinetic properties (/%) for the individual 
binary boundaries, including their anisotropy (see the 
formulation in Sections IV.C and IV.D). 



FIG. 1. (Color online) Time evolution in an TV = 5 
symmetric Cahn-Hilliard problem with starting conditions 
ui — U 2 — Us = U 4 — = 0.2: Different colors stand 

for different phases (where u% > 0.99). Snapshots taken at 
t m (1, 3,10, 30,100, 300) X 25, 000 dimensionless time are dis¬ 
played. (Time increases from left to right, and from top to 
bottom.) 


B. Multiphase separation 

For practical reasons, the time evolution of multi- 
domain systems is investigated in two dimensions using 
conserved dynamics (TV — component Cahn-Hilliard sys¬ 
tems), for which the dimensionless equations of motion 
read as 


Since the dynamic generation of spurious phases is ex¬ 
cluded by the modified Bollada-Jimack-Mullis mobility 
matrix, the trivial extensions of the equilibrium planar 
solutions represent equilibrium, and the stationary solu¬ 
tions of the dynamic equations coincide with the equi¬ 
librium solutions, hence the mobility matrix is regarded 
as ’not dangerous’, when considering planar interfaces. 
Nevertheless, it still remains unclear, whether the same 
applies for all equilibrium solutions, such as the trivial 
extensions of equilibrium trijunctions, etc.. Therefore, 
next we investigate the time evolution of multi-domain 
systems in this respect. 



K \j 


v, £ 



Our choice is motivated by the fact that simulations 
with conserved dynamics are much less forgiving to 
possible numerical errors than those with non-conserved 
dynamics. 


First, we consider a symmetric system, and demon¬ 
strate that our free energy functional prescribes the hier¬ 
archy Fbulk ^ -Winter face F tr i junction <C Fquadruple ^ 
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... for the equilibrium solutions, therefore, a 2- 
dimensional, multi-domain system displays a bulk- 
interface (-trijunction) topology. Accordingly, the ad¬ 
ditional phases near any bulk / interface / trijunction 
vanish with time, since these states are energetically 
not preferred in the vicinity of the equilibrium solu¬ 
tions, to which the system converges. This happens in¬ 
dependently of the particular choice of the mobility ma¬ 
trix ; the only requirement is that the mobility matrix 
has to be symmetric and positive semidefinite. Since 
12 iLi {5F/5ui) = 0 applies for the symmetric model, the 
Lagrangian mobility matrix Kij = 1 yields: 

dUi = 2 SF 

dt 5ui 

a fairly simple system of dynamic equations. Since 
(SF/Sui ) Ui= o = 0, this dynamics satisfies the condition 
of ”no spurious phase generation”. Indeed, we demon¬ 
strate that if a phase is not present in the beginning 
of the calculation, it never appears, even if the other 
phases are not in equilibrium. 

Next, we perform simulations in an N = 4 asym¬ 
metric system with two types of mobility matrices (the 
Lagrangian and the Bollada-Jimack-Mullis type) to 
demonstrate spurious phase generation. We show that 
the spurious phase appears with the Lagrangian mobil¬ 
ity, whereas in the case of the Bollada-Jimack-Mullis 
mobility matrix the zero-amplitude initially prescribed 
for one of the field is retained throughout the simulation, 
even if the other fields are not in equilibrium. Finally, a 
similar behavior is demonstrated for anisotropic systems. 

If not stated otherwise, the dynamic equations were 
solved numerically, using a pseudo-spectral semi-implicit 
method based on operator splitting (see Appendix C), 
on a rectangular grid of size 512 x 512, applying pe¬ 
riodic boundary conditions at the perimeters. The di¬ 
mensionless time and spatial steps were At = 1, and 
Ax = 1. The computations were performed in double 
precision on GTX Titan GPU cards. As starting condi¬ 
tion, we prescribe a spatially homogeneous state contain¬ 
ing equal quantities of all the components, supplemented 
by a small amplitude of flux noise. For this composition, 
equilibrium requires the coexistence of N pure phases. 
Accordingly, the transition process requires the forma¬ 
tion and coarsening of these phases. We present illustra¬ 
tions for symmetric ( N = 5), asymmetric (TV = 4), and 
anisotropic (TV = 4) cases. 



FIG. 2. (Color online) TV = 5 symmetric Cahn-Hilliard 
problem with starting conditions u\ m U 2 — 0.5 and U 3 = 
U 4 — U 5 = 0. (a), (b): Snapshots of the phase fields u\ and 
U 2 at t = 2.5 5 dimensionless time; (c) in the left half panel 
11 ^ 31 + \ua\ + 11 x 51 is shown, whereas the right half panel shows 
m — 1. Note the absence of spurious phases [criterion 
(ii)] and that the sum of the phase fields is 1 [criterion (i)]. 
(d) Multiphase-field map at the same instant. 


0.2,0.2}, which was perturbed by a small amplitude 
of initial noise to induce phase separation. Note the 
coarsening of the various types of grains with time. To 
test whether spurious phase generation takes place, we 
have performed two TV = 5 simulations with initial 
conditions: ( 1 ) u(r, 0) = {0.5,0.5, 0, 0, 0} and for (2) 
u(r, 0) = {0.25, 0.25, 0.25, 0.25, 0}. The individual phase- 
field maps are shown for t = 250, 000 dimensionless time 
in Figs. 2 and 3. Apparently, Uj = 1 is satisfied with a 
high accuracy [criterion (i), see Figs. 2(c) and 3(f)]. Fur¬ 
thermore, the zero amplitude fields retained accurately 
their zero amplitude status throughout the simulations, 
i.e., no third phase generation took place at the phase 
boundaries, indicating that, in agreement with the expec¬ 
tations, criterion (ii) is also fully satisfied. While in the 
effectively two-component case (Fig. 2), we observe the 
usual binary phase separation pattern, the patterns ap¬ 
pearing in Figs. 1 and 3 are significantly different: multi¬ 
grain networks appear that are dominated exclusively by 
trijunctions and binary boundaries; higher-order junc¬ 
tions have not been observed at all. 


1 . Symmetric case 

Snapshots illustrating the time evolution of an TV = 5 
symmetric Cahn-Hilliard problem (where e?- = 1 and 
Wij = 1 for all phase pairs, = 1 to TV) are dis¬ 
played in Fig. 1. The simulation started from a spatially 
homogeneous initial condition u(r, 0) = {0.2,0.2,0.2, 


2. Asymmetric case 

Here, the parameters efj and are different for the 
individual binary interfaces (for the matrices used in the 
present work see Ref.®). The simulations were per¬ 
formed for an TV = 4 Cahn-Hilliard model. u(r, 0) = 
{1/3,1/3,1/3,0} has been chosen as the initial condi- 
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FIG. 3. (Color online) N — 5 symmetric Calm-Hilliard 
problem with starting conditions u\ — 1x2 = 1x3 = 1x4 — 0.25 
and 1 x 5 — 0. (a) - (b): Snapshots of the phase fields ixi, 1x2, 1x3 
and 1 x 4 at t — 10 6 dimensionless time; (e) in the left half panel 
1 x 5 is shown, whereas the right half panel shows X)x=i u i ~ 1 - 
Note the absence of spurious phases [criterion (ii)] and that 
the sum of the phase fields is 1 [criterion (i)]. (f) Multiphase- 
field map at the same instant. 


tion, with a small amplitude of initial noise to induce 
phase separation. First a Lagrangian mobility matrix 
has been used. The phase-field maps corresponding to 
t = 10 6 dimensionless time are displayed in the left col¬ 
umn of Fig. 4. Apparently, the solution is less satisfac¬ 
tory than the results for the symmetric case: substantial 
deviation from 1 x 4 ( 1 *, t) = 0 is observed [Fig. 4(g)]. How¬ 
ever, this spurious phase generation disappears entirely, 
if the Bollada-Jimack-Mullis mobility matrix is used [see 
Fig. 4(h)], as expected. Again, the multiphase domain 
structure is dominated by trijunctions and binary bound¬ 
aries at all times. It appears though that the structure 
obtained with the Bollada-Jimack-Mullis type mobility 
matrix contains chains of alternating u\ and 1 x 4 ” bub¬ 
bles” [see Fig. 4(j)], a feature that can be associated with 
the asymmetry of the kinetic coefficients (/%) applied in 
this simulation. 



FIG. 4. (Color online) N — 4 asymmetric Cahn-Hilliard 
problem with initial condition u = {1/3,1/3,1/3, 0}. The 
results on the left were obtained with a Lagrangian mobility 
matrix, whereas those on the right were obtained using the 
Bollada-Jimack-Mullis type mobility matrix. The upper four 
rows [panels (a) - (h)] show the maps for ixi,iX2,iX3 and 1x4, 
respectively, whereas in the lowermost row [panels (i) and (j)] 
phase distribution maps are displayed. Results corresponding 
to t — 10 6 dimensionless time are presented. Note that using 
the Lagrangian mobility matrix spurious phase generation in 
the vicinity of the phase boundaries could not be avoided for 
1 x 4 [see panel (g)], whereas the Bollada-Jimack-Mullis mobility 
matrix suppresses the formation of spurious phases entirely 
[see panel (h)]. 
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FIG. 5. (Color online) N = 4 asymmetric Calm-Hilliard 
problem with anisotropy: (a) Time evolution towards the 
equilibrium shapes of u\ and U 2 embedded into U 3 at t = 
3.125 X 10 6 dimensionless time, (b)-(e): Snapshots of the 
phase fields ui,U 2 ,U 3 , and 1 x 4 , respectively, taken at t — 
2 x 10 5 dimensionless time. Note the lack of third phase gen¬ 
eration at the phase boundaries, and that ua remains absent 
even at the phase-boundaries. 


3. Anisotropic case 

Next, we investigated an asymmetric N = 4 Cahn- 
Hilliard theor)@, however, with anisotropic interface free 
energies and a Bollada-Jimack-Mullis type mobility ma¬ 
trix. The dynamic equations were solved on a rectan¬ 
gular grid of size 1024 x 512. Time and spatial steps of 
At = 0.25 and Ax = 0.5 were used. The starting con¬ 
ditions were as follows: Two circles filled with - 1 x 1 = 1 
(left) and u 2 = 1 (right) were placed besides each other, 
while a zero value was assigned to these fields outside the 
circles. In contrast, U 3 = 1 was prescribed in the back¬ 
ground, and us = 0 inside the circles, whereas U4 = 0 
was assigned to the whole simulation domain (i.e., the 
fourth field was missing everywhere). All interfaces were 
assumed isotropic, except for the 1-3 interface, for which 
an anisotropy of <213 = 0.1 was prescribed that is larger 
than the critical anisotropy a c = l/(2 fe — 1) = 1/15 for 
fourfold (k = 4) symmetr ^ 49 * 5Q l With elapsing time, 
the circle on the left evolved into a square-like object of 
curved sides, and four pointed corners (see Fig. 5), dis¬ 
playing missing orientations (following from ais > a c ), as 
expected on the basis of the prescribed anisotropy func¬ 
tion. Apparently, as found for the central finite differenc¬ 
ing scheme^, the spectral discretization regularized the 
high anisotropy problem: the predicted numerical shape 
is very close to the analytical solution corresponding to 
this anisotropy. Remarkably, no spurious phase appear¬ 
ance was observed at the phase-boundaries, and = 0 
has been satisfied throughout the simulation. We have 
obtained similar results using finite difference discretiza¬ 
tion. 


C. Grain coarsening 

In this subsection, we apply the XMPF model for grain 
coarsening in a two-dimensional (2D) polycrystalline sys¬ 
tem that contains a large number of differently oriented 
crystal grains that have equal free energy, therefore, the 
time evolution of the system is driven by the grain bound¬ 
ary energy. For the sake of simplicity, we distinguish only 
30 orientations represented by N = 30 fields. The respec¬ 
tive non-conservative equations of motions read as: 



which have been solved numerically on a rectangular grid, 
using a finite difference discretization and Euler forward 
time-stepping. As starting condition, the simulation box 
was covered by a large number of small random field 
patches arranged on a square lattice (512 x 512 and 
256 x 256, respectively, for the larger and smaller size 
simulations), mimicking athermal nucleation on a fine 
grid. [We note that, during time evolution, the initial 
condition is fast forgotten: for example, after a transient 
period, very similar results were obtained with a uniform 
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FIG. 6. (Color online) Snapshot of grain/phase-map taken 
at dimensionless rime, t = 10 4 , for a simulation relying on 
misorientation dependent grain boundary energy obeying the 
Read-Shockley relationship 51 . The simulation was performed 
on a 4096 2 rectangular grid. About 4000 grains can be dis¬ 
tinguished at this stage. 




r/<r> 


1 /N starting, while adding a small pixelwise Gaussian 
noise (a spatially random initiation).] Two cases were in¬ 
vestigated: (a) with an isotropic grain boundary energy 
(on a 8192 2 grid), and (b) the misorientation dependence 
of the grain boundary energy follows the Read-Shockley 
relationship 51 ! (on a 4096 2 grid). The corresponding mo¬ 
bilities were K>ij = 1, and Kij = \ui/(\ — Uj) \ \uj/(l — Uj)\, 
respectively. 

A typical grain map displaying the result of the simu¬ 
lation for case (b) at a dimensionless time t = 10 4 , when 
~ 4000 grains exist, is shown in Fig. 6. Similar grain 
maps have been obtained for the other case, except that 
there most of the trijunctions are close to symmetric, 
displaying angles ~ 120°. 

As observed in other MPF model^^J, in the 
experiments^, and predicted by theor}®^, after a tran¬ 
sient period a limiting grain size distribution (LGSD) is 
established, which in the case of experiments on metallic 
film can be accurately fitted^ by the lognormal distri¬ 
bution proposed by Felthamf^. The models of Mullins 
and Hillert predict significantly different LGSDs (Fig. 7). 
The LGSD predicted by the XMPF model approximates 
the experimenta l resu lts somewhat better than the previ¬ 
ous MPF models 27 ^ [see Figs. 7(b) and 7(c)], and prac¬ 
tically coincides with the re sults from the mutli order pa¬ 
rameter approaches 11 ^ yet the agreement is not par¬ 
ticularly good with the experiments at small grain sizes. 
Apparently, in the experiments the small grains disap- 


FIG. 7. (Color online) Limiting grain size distributions 
(LGSD): (a) for the isotropic simulation on a 8192 2 grid. For 
comparison, experimental result^ 53 for metallic films (solid 
line: lognormal distribution fitted to experimental dateFp, 
and predictions by the theories of Mulling and HillertP^ 
(dashed and dotted lines, respectively) are also shown, (b) 
Comparison of the LGSD from XMPF (histogram) with pre¬ 
dictions from previous MPFs by Kim et air 7 ' (triangles) and 
Schaffnit et a/“ (circles). The solid line indicates a lognor¬ 
mal fit to the experimental results 53 } (c) The same as panel 
(b), except that LGSD from an asymmetric XMPF model 
(the simulation shown in Fig. 6) is displayed (histogram), in 
which the grain boundary energy follows the Read-Shockley 
relationship 5 ^. Note that some improvement relative to the 
previous MPF models has been achieved, but the population 
of the small grains is larger than desirable. 


pear faster than in the XMPF simulations. In the investi¬ 
gated cases, the time dependence of the average grain size 
[(r) = A(t—to) q , where A is a constant and to the freezing 
time] is described by an exponent q = 0.5(1 ± 0.05), in¬ 
dicating an essentially diffusion controlled grain growth. 

We mention in this respect that a simple dynami¬ 
cal density functional theory, the Phase-Field Crystal 
approacfPH, which incorporates a broad range of phys¬ 
ical phenomena (elasticity, dislocation dynamics, grain 
rotation, etc.), reproduces the experimental LGSD fairly 
welP^. Unfortunately, in the PFC studies, as in the case 
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of experiments, the effect of different physical phenomena 
on the LGSD cannot be easily separated. It is expected, 
however, that the comparison of different models may 
contribute to the identification of the governing phenom¬ 
ena. Along these lines, the present study determined 
the LGSD the physically consistent XMPF model pre¬ 
dicts. Apparently, further efforts are needed to improve 
the agreement between MPF models and experiments. 
Work is, underwa^ 9 ! to evaluate LG SD fro m phase-field 
models relying on orientation field (s)l^® 0 ^ in describing 
different crystallographic orientations. 


VI. COMPARISON WITH OTHER MODELS 

Having presented the essential properties of the 
XMPF model, it is desirable to compare it with other 
models from the following viewpoints: 

A. For a few of the most important multiphase-field 
models, we investigate whether the trivial N = 3 exten¬ 
sion of the equilibrium binary solution is a stationary 
solution of the N = 3 dynamic problem [part of criterion 
(vi), which in addition requires that the trivial extension 
be a solution of the N = 3 Euler-Lagrange equation 
too]. 

B. We explore, furthermore, for the best behaving 
models identified in sub-section A. whether the free en¬ 
ergy decreases indeed monotonically with time [criterion 

(v)]. 


A. Planar interfaces: 
comparing the MPF models 

Here, we investigate for several MPF models, whether 
a trivial N = 3 extension of the equilibrium interface 
of the two-phase problem (obtained by adding ^3 = 0 
to the two-phase equilibrium solution) behaves like a lo¬ 
cal free energy minimum of the multiphase-field model: 
starting from the extended solution, we explore whether 
the N = 3 equations of motion keep the solution equal 
to the starting condition, or drive it away. For this test, 
we adopt non-conservative dynamics 



The initial condition is a liquid-solid-liquid slab, with 
periodic boundary condition at the two ends, while 
employing the respective analytic solutions in the inter¬ 
face regions, accompanied with Us(z) = 0 throughout 
the computation box. (See Figs. 8(a) and 9(a) for the 
initial conditions used for the models that have binary 
equilibrium solutions of the forms u(z) = [1 + sin(z)]/2 




FIG. 8. (Color online) u(z) = [1 + sin(z)]/2 models, (a) 
Initial condition (top panel), (b)-(g) final states after 10 6 time 
steps. From second to bottom row, respectively: the model 
of Nestler and Wheeler p — l 17 19 20 ^, the model of Steinbach 
and Pezzolgp^, and the model of Steinbach and Pezzola with 
non-variational dynamics 29 . (Left column: final states, right 
column: difference of final state and initial condition.) 


and u(z) = [1 + tanh(z/2)]/2, respectively). 


The one-dimensional dynamic equations were solved 
numerically, using finite difference method, while em¬ 
ploying dimensionless time and spatial steps, h and At, 
as specified below. 
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FIG. 9. (Color online) u(z) — [1 + tanh(z/2)]/2 mod¬ 
els. (a) Initial condition (top panel), (b)-(g) final states after 
2.5 x 10 5 time steps. From second to bottom row, respectively: 
the model of Steinbach et (coincides with the model of 
Nestler and Wheeler for p = 2 16 ^), the model of Steinbach 
et al. with non-variational dynamics 1 ^, and the XMPF model 
proposed in this work. (Left column: final states, right col¬ 
umn: difference of final state and initial condition.) 


1. Models with sinusoidal equilibrium profile 

The long-time solutions of the dynamic equations (10 6 
time steps, beyond which no further changes were per¬ 
ceptible) are shown in Fig. 8 for the Nestler-Wheeler 
p = 1 modeff^^Sl[Figs. 8(b) and 8(c)], for the 

Steinbach-Pezzola model 13 ' [Figs. 8(d) and 8(e)], and 
for the Steinbach-Pezzola model with non-variational 
dynamicd^ [Figs. 8(f) and 8(g)] (h = 0.1 and At = 


0.0025 were used.) The long-time interfacial field pro¬ 
files are shown on the left [Figs. 8(b), 8(d), and 8(f)], 
together with their difference relative to the initial con¬ 
ditions on the right [Figs. 8(c), 8(e), and 8(g)]. While in 
the first and second cases, third-phase generation can be 
seen at the interface, the application of non-variational 
dynamics in the Pezzola-Steinbach model suppressed this 
phenomenon (us ~ 0 was retained). 


2. Models with hyperbolic tangential equilibrium profile 

The long-time solutions (2.5 x 10 5 time steps, beyond 
which no perceptible changes were seen) of the EOMs 
are shown in Fig. 9 for the following models: Figs. 9(b) 
and 9(c) - the model of Steinbach et al.^\ Figs. 9(d) and 
9(e) - the model of Steinbach et al. ^ with non-variational 
dynamics; Figs. 9(f) and 9(g) - the model proposed in 
the present paper, (h = 0.25 and At = 0.01.) The 
long-time interfacial field profiles are shown on the left 
[Figs. 9(b), 9(d), and 9(f)], together with their difference 
relative to the initial conditions on the right [Figs. 9(c), 
9(e), and 9(g)]. While the model of Steinbach et al 12 
leads to third-phase generation, the other two approaches 
are free of this problem. Remarkably, the predictions 
from the latter two models fall very close to each other. 
Yet, in the model of Steinbach et al the trivial three- 
phase extension of the binary equilibrium solution is not 
a solution of the three-phase Euler-Lagrange equation 
(see Section III.C). In other words: although the same 
solution is a stationary solution of the non-variational 
EOM, stabilized by the non-variational dynamics, it is 
not a free energy minimum of the three-phase problem. 

While in this test, the results of the model of Stein¬ 
bach et al. (with non-variational dynamics) are practi¬ 
cally indistinguishable from those of the XMPF model 
proposed in this work, under other conditions significant 
differences can be seen. 


B. Time dependence of the free energy 

In this test, we investigate the time evolution of the 
system in the non-variational model of Steinbach et al 
and in the XMPF model presented in this work. A sym¬ 
metric N = 4 Cahn-Hilliard model and Lagrangian mo¬ 
bility matrix have been chosen for the demonstration. 

The results are summarized in Fig. 10, which shows 
the map of one of the fields (the others are qualitatively 
similar) at dimensionless times t = 2000 and 20000, com¬ 
puted on a grid 512 x 512. While in the XMPF model 
proposed here, the free energy decreases monotonically 
with time as expected, in the model of Steinbach et al. 
(when relying on a non-variational formalism), the free 
energy increases initially, reaching then a maximum at 
about t = 3000, followed by a slow decrease beyond. 
This behavior is presumably a consequence of the ap¬ 
plied ’binary’ approximation, in which various terms of 
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FIG. 10. (Color online) Comparison of the time evolution 
of phase transition (a), (c) in the non-variational model of 
Steinbach et (left), and (b), (d) in the present model 

(right). Snapshots of one of the fields are displayed. The time 
dependence of the free energy is shown in panel (e): dashed 
line - model of Steinbach et a/.; solid line - the XMPF model. 
The symbols correspond to snapshots shown in panels (a)-(d). 


the variational equations of motion are omitted: once the 
free energy functional (Lyapunov functional) is defined, 
the variational dynamics ensures a monotonic reduction 
of the free energy with time. Any deviation from this 
approach raises the possibility of a non-monotonic time 
evolution of the free energy. 


VII. SUMMARY 


In this work, we formulated a physically consistent 
multiphase-field theory for describing interface driven 
multi-domain processes. First, we identified a set of crite¬ 
ria, a physically consistent multiphase-field approach has 
to satisfy. These are: (i) the sum of the fields is 1 every¬ 
where; (ii) the physical results should be invariant to ex¬ 
changing pairs of field indices, i j\ (iii) a trivial multi¬ 
phase extension of the equilibrium binary solution should 
represent an equilibrium solution of the multiphase prob¬ 
lem, which in turn should be a stationary solution of the 
dynamic equations towards which the time dependent so¬ 
lutions evolve; (iv) variational dynamic equations shall 


be used to ensure non-negative entropy production; (v) 
reduction/extension of the TV-held theory to N — 1 or 
N + 1 fields should be straightforward, and happen con¬ 
sistently within the formalism; (vi) there should be no 
spurious third phase appearance at the equilibrium bi¬ 
nary boundaries, and once a field is not present, it should 
not appear at any time in the dynamic equations; finally, 
(vii) freedom to choose the interfacial and kinetic prop¬ 
erties for individual phase pairs. 

Next, considering these requirements, we have re¬ 
viewed a range of the existing multiphase field models, 
and identified their advantageous and less advantageous 
features. 

Combining the advantageous features of the ear¬ 
lier multiphase-field models, we have constructed a 
multiphase-field approach (termed the XMPF model) 
obeying all criteria defined above. In addition, we 
performed illustrative simulations for N = 4 and 5 
multiphase-field models that rely on conserved dynam¬ 
ics, describing thus multiphase separation problems (TV- 
component Cahn-Hilliard problems). Symmetric (identi¬ 
cal interface properties), asymmetric (pairwise different 
interface properties), and anisotropic (orientation depen¬ 
dent interfacial properties) cases were addressed, and it 
has been shown that using a suitable mobility matrix 
(Bollada-Jimack-Mullis type), the XMPF model avoids 
dynamic spurious phase generation. We have performed 
further illustrative simulations for grain coarsening in 
polycrystalline systems using an N = 30 XMPF model 
relying on non-conserved dynamics. While the predicted 
limiting grain size distribution is closer to the experimen¬ 
tal results than those from the previous MPF models, 
further works is needed to improve the agreement. 

The present work opens up the way towards physically 
consistent computations for microstructure evolution in 
multiphase / multigrain / multicomponent structures, 
and shall serve as a basis for developing a physically con¬ 
sistent quantitative multiphase-field approach that might 
be combined with melt flow and elasticity, and extended 
to fast processes along the lines described in RefsP^®, 
leading towards developing improved tools for knowledge 
based materials design. 

Work is underway to incorporate a phase-dependent 
thermodynamic driving force (a multiphase analogy of 
the ’tilting function’ in RefP^) into the XMPF model, 
which will be presented in a separate paper. 

We note in this respect that the inclusion of thermo¬ 
dynamic driving force via a tilting function has no effect 
on the present results concerning the two- and multi¬ 
phase equilibria. The existence of equilibrium two-phase 
planar interfaces in the multiphase problem is a basic 
requirement, which needs to be satisfied by a physically 
consistent model. 
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APPENDIX Al: INVARIANCE OF RESULTS TO 
EXCHANGING PAIRS OF FIELD INDICES, i ^ j 

The general dynamic equations of a multiphase model 
read as: 



where there are N(N — 1) mobilities (/%). The princi¬ 
ple of formal indistinguishability of the variables means 
that the variables are not ”labeled”, i.e. none of them 
is distinguished formally on the basis of its index. This 
is true if the dynamic equations are invariant for the re¬ 
labeling of the variables, i.e. re-labeling of the variables 
on the level of the free energy functional results in the 
same as re-labeling the variables in the dynamic equa¬ 
tions. This criterion is satisfied by symmetric mobility 
matrices, namely, 


Hvij ^ji • 


Proof. The dynamic equation for uj reads as 


duj 

dt 


^ Kjk 
k^J 


f SF 
\Suj 



(38) 


The variables can be re-labeled by using the variable 
transformation v & := u & for k ^ I, J, Vi := uj , and 
vj := ui. Using this in Eq. (38) yields then 


dvj 

~dt 



SF _ SF\ 
Svi Svj J 


where the chain rule for the functional derivative has also 
been used (see Appendix A2). Furthermore, re-labeling 
variables in the free energy functional first (F[ u] —>■ 
F[v]), then deriving the dynamic equations simply re¬ 
sults in 



Comparing the two equations yields then ku — kj/. 


In order to illustrate the ”no labeling” condition in 
practice, we choose a typical example of labeling the vari¬ 
ables. Some authors eliminate of one of the variables even 


at the level of the free energy functional, i.e. they intro¬ 
duce the independent variables Vi := tq for i = 1... N— 1, 
thus resulting in un = 1 — v^. Then, the following 

dynamic equations are used: 

dvi SF 

T% dt Svi 

These can be written in terms of the old variables as: 


dui _ 1 / SF SF \ 

dt Ti \Sui Sun J 

for i = 1... N — 1, and 

SF \ 

Su N J 



(39) 


(40) 


It is straightforward to see, that Eqs. 
scribe the following mobility matrix: 


(39) and (40) pre- 


Tii — 


and LiN = —r- 


for i = 1... N — 1, while the last row reads as 


N—l 


Lni = —T- 


-1 


and L nn = ^ r { 1 


i= 1 


where the form —dui/dt = Lij{SF/S u j ) is used. It 

is trivial that the elements of L sum up to 0 in each row 
and column, but the matrix is not symmetric! It means 
that the concept of eliminating a variable on the level of 
the free energy functional labels the variables, i.e. the 
eliminated variable is formally distinguished. Indeed, 
exchanging variables I and V, deriving the dynamic 
equations, then exchanging I and N back result in a 
mobility matrix similar to the one described by Eqs. 
(39) and (40), however, the N th and the I th rows are 


exchanged. On the one hand, it means that the for¬ 
mal variable exchange corresponds to the elimination of 
phase I instead of phase N. On the other hand, since the 
resulting mobility matrix is not identical to the original 
one, the eliminated variable is always labeled, therefore, 
this concept does not satisfy the condition of no labeling. 


APPENDIX A2: CHAIN RULE FOR 
FUNCTIONAL DIFFERENTIATION 

Mathematically speaking, the solution of the Euler- 
Lagrange equations is invariant to the variable transfor¬ 
mation Q = T[q], if the transformation T[.] is unam¬ 
biguous, i.e., if the inverse transform T -1 [.] also exists. 

Proof. The Euler-Lagrange equations for the new vari¬ 
ables read as: 


SF_dI dl 

SQi dQi d\/Qi 


( 41 ) 
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where I denotes the full integrand of Eq. <©• The terms 
on the right-hand side can be expanded as follows: 


dl ^ dl dqj dl d\/qj 

dQi = ^d^jdQ~i + 9Vqj~dQT 


V 


dVQi 


= V 


E 


dl dqj 
dqj dVQi 


dl dVgj 
dVqj dVQi 


.(43) 


Since qj = T- (Q), dqj/d\7Qi = 0. In addition, formally 
V <7j = Ei §o~ V Qi’ therefore, S* = lY Using these 


together with Eqs. (42) and (43) in Eq. (41) yields 


APPENDIX B: NUMERICAL METHOD 

The dynamic equations were solved numerically on a 
periodic, two-dimensional domain by using an operator- 
splitting based, quasi-spectral, semi-implicit time step¬ 
ping scheme as follows. The dynamic equations can be 
re-written in the form 

^ = f(u,Vu), (46) 

where f(u, Vu) is the general, non-linear right-hand side. 
During time stepping f(u, Vu) is calculated at time point 
£, while dui/dt is approximated as 


iL = W—-V—i) 

8Qi “ \dqj d\7qj) dQi 

91 ( (Ngj _ y 9<lj_ 

A dVqj V dQi do, 


Finally, Vqj = Efc dOfc^Qk -» r ~qq~ — L~,k OQ^dQ. 
and V§q~ = Efe 3 Q VQfc, therefore, the second sum 
on the right hand side of Eq. (44) vanishes. The final 
result then reads as: 




dQi 




JLL 


VQ/c 


SF dqj 

$Qi . dQi 


(45) 


i.e. the chain rule of differentiation also applies for the 
functional derivative. Let now q* (r) denote the solution 
of SF/S q = 0. Apparently, the right hand side of Eq. 
( [45] ) vanishes at q*(r). Formally q*(r) = T _1 [Q*(r)], 
indicating that Q*(r) = T[q*(r)], i.e. the solution in Q 
is just the transformation of the solution in q. In other 
words, the solution of the Euler-Lagrange equations is 
invariant to the choice of the generalized variables. 


duj 

dt 


■■r A< - < 

At 


(47) 


(44) Next, we add a suitably chosen linear term s[ui] = 


• 1 (—1 fsiW^Uj (where s* > 0) to both sides of Eq. 
( [46] ) . We consider this term at t + At at the left-hand 
side, but at t on the right-hand side of the equation. This 
concept, together with Eq. (47) results in the following, 


explicit time stepping scheme in the spectrum: 


4+At 


At 


(k) = ^ (k) + TT^(k)^^ [u(r) ’ Vu(r)]} 


where s^(k) = JfJLi s^(k 2 ) J *, and F{.} stands for the 
Fourier transform. The splitting constants must 


(48) 


be chosen so that Eq. (48) to be stable 


It is important to note that our numerical scheme is 
unbounded , which means that the spatial solution rq(r, t) 
can go under 0 or above 1 because of the numerical errors. 
The construction of the free energy functional and the 
modified Bollada-Jimack-Mullis mobility matrix, how¬ 
ever, ensure that the system converges to equilibrium. 
This means that no artificial modification of the solution 
is needed after a time step, which could lead to instabil¬ 
ities in the spectral method. 
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